{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data poisoning attack\n",
    "\n",
    "In this notebook, we use a convex optimization layer to perform a *data poisoning attack*; i.e., we show how to perturb the data used to train a logistic regression classifier so as to maximally increase the test loss. This example is also presented in section 6.1 of the paper *Differentiable convex optimization layers*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cp\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "\n",
    "\n",
    "from cvxpylayers.tensorflow.cvxpylayer import CvxpyLayer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are given training data $(x_i, y_i)_{i=1}^{N}$,\n",
    "where $x_i\\in\\mathbf{R}^n$ are feature vectors and $y_i\\in\\{0,1\\}$ are the labels.\n",
    "Suppose we fit a model for this classification problem by solving\n",
    "\\begin{equation}\n",
    "\\begin{array}{ll}\n",
    "\\mbox{minimize} & \\frac{1}{N}\\sum_{i=1}^N \\ell(\\theta; x_i, y_i) + r(\\theta),\n",
    "\\end{array}\n",
    "\\label{eq:trainlinear}\n",
    "\\end{equation}\n",
    "where the loss function $\\ell(\\theta; x_i, y_i)$ is convex in $\\theta \\in \\mathbf{R}^n$ and $r(\\theta)$ is a convex\n",
    "regularizer. We hope that the test loss $\\mathcal{L}^{\\mathrm{test}}(\\theta) =\n",
    "\\frac{1}{M}\\sum_{i=1}^M \\ell(\\theta; \\tilde x_i, \\tilde y_i)$ is small, where\n",
    "$(\\tilde x_i, \\tilde y_i)_{i=1}^{M}$ is our test set. In this  example, we use the logistic loss\n",
    "\n",
    "\\begin{equation}\n",
    "\\ell(\\theta; x_i, y_i) = \\log(1 + \\exp(\\beta^Tx_i + b)) - y_i(\\beta^Tx_i + b)\n",
    "\\end{equation}\n",
    "\n",
    "with elastic net regularizaiton\n",
    "\n",
    "\\begin{equation}\n",
    "r(\\theta) = 0.1\\|\\beta\\|_1 + 0.1\\|\\beta\\|_2^2.\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.datasets import make_blobs\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "\n",
    "tf.random.set_seed(0)\n",
    "np.random.seed(0)\n",
    "n = 2\n",
    "N = 60\n",
    "X, y = make_blobs(N, n, centers=np.array([[2, 2], [-2, -2]]), cluster_std=3)\n",
    "Xtrain, Xtest, ytrain, ytest = train_test_split(X, y, test_size=.5)\n",
    "\n",
    "Xtrain, Xtest, ytrain, ytest = map(\n",
    "    tf.constant, [Xtrain, Xtest, ytrain, ytest])\n",
    "m = Xtrain.shape[0]\n",
    "\n",
    "lambda1_tf = tf.constant([[0.1]], dtype=tf.float64)\n",
    "lambda2_tf = tf.constant([[0.1]], dtype=tf.float64)\n",
    "\n",
    "a = cp.Variable((n, 1))\n",
    "b = cp.Variable((1, 1))\n",
    "lambda1 = cp.Parameter((1, 1), nonneg=True)\n",
    "lambda2 = cp.Parameter((1, 1), nonneg=True)\n",
    "X = cp.Parameter((m, n))\n",
    "Y = ytrain.numpy()[:, np.newaxis]\n",
    "\n",
    "log_likelihood = (1. / m) * cp.sum(\n",
    "    cp.multiply(Y, X @ a + b) -\n",
    "    cp.log_sum_exp(cp.hstack([np.zeros((m, 1)), X @ a + b]).T, axis=0,\n",
    "                   keepdims=True).T\n",
    ")\n",
    "regularization = - lambda1 * cp.norm(a, 1) - lambda2 * cp.sum_squares(a)\n",
    "prob = cp.Problem(cp.Maximize(log_likelihood + regularization))\n",
    "fit_logreg = CvxpyLayer(prob, [X, lambda1, lambda2], [a, b])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume that our training data is subject to a data poisoning attack,\n",
    "before it is supplied to us. The adversary has full knowledge of our modeling\n",
    "choice, meaning that they know the form of the optimization problem above, and seeks\n",
    "to perturb the data to maximally increase our loss on the test\n",
    "set, to which they also have access. The adversary is permitted to apply an\n",
    "additive perturbation $\\delta_i \\in \\mathbf{R}^n$ to each of the training points $x_i$,\n",
    "with the perturbations satisfying $\\|\\delta_i\\|_\\infty \\leq 0.01$.\n",
    "\n",
    "Let $\\theta^\\star$ be optimal.\n",
    "The gradient of\n",
    "the test loss with respect to a training data point, $\\nabla_{x_i}\n",
    "\\mathcal{L}^{\\mathrm{test}}(\\theta^\\star)$, gives the direction\n",
    "in which the point should be moved to achieve the greatest\n",
    "increase in test loss. Hence, one reasonable adversarial policy is to set $x_i\n",
    ":= x_i +\n",
    ".01\\mathrm{sign}(\\nabla_{x_i}\\mathcal{L}^{\\mathrm{test}}(\\theta^\\star))$. The\n",
    "quantity $0.01\\sum_{i=1}^N \\|\\nabla_{x_i}\n",
    "\\mathcal{L}^{\\mathrm{test}}(\\theta^\\star)\\|_1$ is the predicted increase in\n",
    "our test loss due to the poisoning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "\n",
    "loss = tf.keras.losses.BinaryCrossentropy()\n",
    "with tf.GradientTape() as tape:\n",
    "    tape.watch(Xtrain)\n",
    "    # Apply the layer\n",
    "    slope, intercept = fit_logreg(Xtrain, lambda1_tf, lambda2_tf)\n",
    "    # 30 is scale factor so visualization is pretty\n",
    "    test_loss = 30 * loss(ytest, Xtest @ slope + intercept)\n",
    "# Compute the gradient of the test loss with respect to the training data\n",
    "Xtrain_grad = tape.gradient(test_loss, Xtrain)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below, we plot the gradient of the test loss with respect to the training data points. The blue and orange points are training data, belonging to different classes. The red line is the hyperplane learned by fitting the the model, while the blue line is the hyperplane that minimizes the test loss. The gradients are visualized as black lines, attached to the data points. Moving the points in the gradient directions torques the learned hyperplane away from the optimal hyperplane for the test set."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hU1dbH8e/OJCEJXar0oICgFCEUsYGAlyZVURBURJqI2BFR0VdRLFcUkas0excQQVCKRlQsgASkg4D0TigphGTW+8dOImogZcqZmazP8+SBTGbOWQnhl5N99l7biAhKKaVCR5jTBSillPIuDXallAoxGuxKKRViNNiVUirEaLArpVSI0WBXSqkQ45VgN8bca4xZa4xZY4z50BgT5Y3jKqWUyj+Pg90YUxm4G4gTkUsAF3CTp8dVSilVMN4aigkHoo0x4UAMsMdLx1VKKZVP4Z4eQER2G2NeBHYAKcACEVnwz+cZYwYBgwCKFi3a5KKLLvL01CFFBP74A44dg2rVoFw5P5788BZISwZxgwmDyBgoc6EfC1BK5cWKFSsOiUiu6WA8bSlgjCkNzABuBBKBT4HPROS9s70mLi5Oli9f7tF5Q9GpU9CzJ3z5JbzxBgwa5IeTbvwKZtyO+9RJ0t0Q6TIQWRR6Toc67f1QgFIqr4wxK0QkLrfneWMopi2wTUQOishpYCbQ0gvHLXSKFIEZM6BjRxg8GKZM8cNJ961m96GTuP7vBJNXpNnH0pJh3+9+OLlSyhe8Eew7gBbGmBhjjAHaAOu9cNxCKSvcO3SwV+xTp/r2fOuTSlFl/AkABjeJtA9GxkDF+r49sVLKZzwOdhH5BfgM+A34PfOYkz09bmEWFQUzZ0L79jBwIEyf7pvz/Pjjj9TrNBgA99MViXCF2WGYynFQq51vTqqU8jmPb54CiMgYYIw3jqWsqCiYNQu6dYM77gBjoH9/7x3/888/p3v37lStWpUd27fB5oV2+KVifRvqYS7vnUwpLzh9+jS7du0iNTXV6VJ8LioqiipVqhAREVGg13sl2JVvREXB559D164wYIAN99tu8/y4r7/+OkOHDqV169Z888039sE67fVmqQpou3btonjx4tSoUQM76huaRITDhw+za9cuYmNjC3QMbSkQ4LLCvW1buP12ePttz4732GOPMXToUPr37/9XqCsVBFJTUylTpkxIhzqAMYYyZcp49JuJBnsQiI6G2bNtuPfvD++8U7Dj3HbbbTz99NOMGTOG6b4auFfKh0I91LN4+nnqUEyQyAr3Ll3scIwx0K9f3l9/1VVX8f333zNlyhTuuOMOn9WplHKeXrEHkaxwv+YauPVWeO+sS8D+rlKlSnz//fd88cUXGupKeSAxMZFJkybl+3UdO3YkMTHRBxXlTIM9yMTEwBdfQOvWNtzff//cz69fvz579+5l6dKlXHfddf4pUqkQdbZgz8jIOOfr5s2bR6lSpXxV1r/oUEwQiomBOXOgc2e45RY7LNOnT87Pfe2116hQoQJ16tTxb5FKhaCHH36YP/74g0aNGhEREUGxYsU4//zzSUhIYN26dXTr1o2dO3eSmprKiBEjGJTZF6RGjRosX76ckydP0qFDB6644gqWLl1K5cqVmT17NtHR0V6tU4M9SJ0Z7v362XDv3fvfz7vqqqv8X5xSvnbPPZCQ4N1jNmoEL798zqeMGzeONWvWkJCQQHx8PJ06dWLNmjXZ0xKnT5/OeeedR0pKCk2bNqVnz56UKVPmb8fYvHkzH374IVOmTKFXr17MmDGDvn37evVT0aGYIFa0KMydC1deCX37wscfO12RUoVLs2bN/jbXfMKECTRs2JAWLVqwc+dONm/e/K/XxMbG0qhRIwCaNGnC9u3bvV6XXrEHuaJFbTfIjh3h5pvtlXuvXk5XpZSP5XJl7S9FixbN/nt8fDyLFi3ip59+IiYmhlatWuU4F71IkSLZf3e5XKSkpHi9Lr1iDwFZ4d6ypR1r/+QTpytSKjQVL16cEydO5PixY8eOUbp0aWJiYtiwYQM///yzn6v7i16xh4hixWDePNsVsk8fCAuD6693uiqlQkuZMmW4/PLLueSSS4iOjqZChQrZH2vfvj2vv/46DRo0oE6dOrRo0cKxOj3eaKMg4po0keUrVvj9vIXBiRM23H/+2Y659+zpdEVKecf69eupW7eu02X4TU6frz832si/tWvhww8hl7mfKv+KF4f586F5c7jpJtv+VylVuDgT7GFhdrygfn07IOx2O1JGqMoK96ZN4cYbbftfpVTh4Uyw16tnA90YmzwNG9pLSw14rylRAr76CuLi7CyZzz93uiKllL84Nyvmhhtg9Wr44AM4fdoOBjdubJuhODDuH4qywr1JE/vlnj3b6YqUUv7g7HRHl8sul1y7Ft59F5KS7JZBTZva+Xsa8B4rWRK+/vqvcP/iC6crUkr5mleC3RhTyhjzmTFmgzFmvTHmsnwdwOWySyfXr4c334QjR+xa+RYtbCppwHskK9wvvdROgZwzx+mKlFK+5K0r9leAr0TkIqAhsL5ARwkPt83GN26EyZNh3z67o/MVV8CiRRrwHsgK90aN7KjX3LlOV6RU8Clo216Al19+meTkZC9XlDOPg90YUwK4CpgGICJpIuJZ4+GICBg4EDZvhv/9D/78E9q1g1at4LvvPC250CpVChYssPeqe/a0o11KqbwrNMEO1AQOAm8aY1YaY6YaY4rm9qI8iYyEIUNgyxaYMMEGfatW0KYN/PCDV05R2GSFe4MG0KOHXa2qlMqbM9v2Pvjgg7zwwgs0bdqUBg0aMGbMGACSkpLo1KkTDRs25JJLLuHjjz9mwoQJ7Nmzh9atW9O6dWuf1+mNlgLhQGNguIj8Yox5BXgYeOzMJxljBgGDAKpVq5a/M0RFwfDhcMcd8MYbMG6cbWnYrh08+SRclr8h/cKudGkb7u3aQffudipkhw5OV6VU3jnUtfdvbXsXLFjAZ599xq+//oqI0KVLF5YsWcLBgwepVKkSX2b+Snzs2DFKlizJSy+9xLfffkvZsmW9W3gOvHHFvgvYJSK/ZL7/GTbo/0ZEJotInIjElStXrmBnio62/6Jbt8KLL8LKlbbzVceOsGxZgT+Bwqh0aVi4EC65xIb7V185XZFSwWXBggUsWLCASy+9lMaNG7NhwwY2b95M/fr1WbRoESNHjuT777+nZMmSfq/N4yt2EdlnjNlpjKkjIhuBNsA6z0s7h5gYuP9+GDwYJk6EF16AZs3guuvsFfyll/r09KEiK9zbtrWzTGfPhv/8x+mqlMpdIHTtFRFGjRrF4MGD//WxFStWMG/ePEaNGsW1117L448/7tfavDUrZjjwvjFmNdAIeMZLxz23YsXg4Ydh2zZ4+mn4/nu7yKlHD7v4SeXqvPPshKO6daFrVztEo5TK2Zlte//zn/8wffp0Tp48CcDu3bs5cOAAe/bsISYmhr59+/LAAw/w22+//eu1vuaVYBeRhMxhlgYi0k1EjnrjuHlWogSMHg3bt8MTT8DixXbqR69edvGTOqescL/oIhvuCxc6XZFSgenMtr0LFy6kT58+XHbZZdSvX5/rr7+eEydO8Pvvv9OsWTMaNWrE2LFjefTRRwEYNGgQHTp08MvNU2fa9sbFyfLly313gqNH4aWX7O9rSUm2zeHjj9vkUmd16JCdcLRpk13E1Lat0xUp9Rdt2xvobXt9rXRpeOopewU/cqRdR3/xxXbX5xz2IFRW2bL2l51ateztisWLna5IKVUQoRnsWcqUgWeftbNo7rsPZsywg8n9+9vH1L/8M9y/+cbpipRS+RXawZ6lfHk7c2brVjsf/qOPoHZtOy/+zz+dri7glCtnw/2CC2zLHg13FSicGDp2gqefZ+EI9iwVK8L48fDHHzB0qO0oWauW/fvOnU5XF1Cywr1mTRvu8fFOV6QKu6ioKA4fPhzy4S4iHD58mKioqAIfIzRvnubVrl3wzDMwdard9GPQIBg1CipVcrqygHHgALRubW9XzJsHV1/tdEWqsDp9+jS7du0iNTXV6VJ8LioqiipVqhAREfG3x/N687RwB3uWP/+0AT99um0hPGSInR9fsaLTlQWE/fvhmms03JVyWuGeFZNf1avbHjSbNsHNN9vVrDVrwgMP2EvWQq5CBTvOXr267d6wZInTFSmlzkWD/UyxsTBtGmzYYLcbGj/ePjZypJ3kXYhlhXu1ajbcv//e6YqUUmejwZ6TCy+Et9+GdevsUswXXrAB/+ijdnenQqpiRfj2W6ha1XaD1M7JSgUmDfZzqVPHbra9Zo29TB071gb8E09Aomd7iQSrihXtlXuVKhruSgUqDfa8qFcPPv7YNhbL6gEfG2tXtx4/7nR1fnf++fbKvVIlG+5LlzpdkVLqTBrs+VG/Pnz2me0Df/XVtv9MbKxd3ZrZ4a2wyAr388+329L+9JPTFSmlsmiwF0SjRnbboWXL7O5NjzxiA/6FF2zTsUKiUiUb7hUr2j7uP//sdEVKKdBg90xcHMydaxOtcWN46CE7TXL8eEhJcbo6v6hc2YZ7hQo23H/5JffXKKV8S4PdG5o3h6+/tncS69e3DccuuABefRUKwSq5rHAvVw6uvRZ+/dXpipQq3DTYvenyy+2OFd99Z5uM3X23nTo5aRKcOuV0dT5VpYoN97JlbbjrFrRKOUeD3Reuusqm3OLFUKMGDBtmm41NngxpaU5X5zNVq9pmYWXK2MlDGu5KOUOD3VeMsQ1Wvv/eDtNUqmQ3365Tx/akOX3a6Qp9ompV+zPtvPNsuAdSSyClCguvBbsxxmWMWWmMmeutY4YEY+zYxE8/2Q5aZcvCgAF2w4933oH0dKcr9Lpq1Wy4ly5tw33FCqcrUqpw8eYV+whgvRePF1qMsat5fv0VZs+G4sXh1lvtln0ffAAZGU5X6FXVq9thmVKl7N6pGu5K+Y9Xgt0YUwXoBEz1xvFCmjHQpYtNupkzoUgR21Gyfn345BNwu52u0GuqV7dX7iVL2iv3335zuiKlCgdvXbG/DDwEnDWVjDGDjDHLjTHLDx486KXTBrGwMOjeHRISbKAbAzfeCA0b2sAPkYCvUcNeuRcvbq/cV650uiKlQp/HwW6M6QwcEJFz/rItIpNFJE5E4sqVK+fpaUNHWJhtEbx6tR2SOX0aeva0C55mz4YQ2Absn+GekOB0RUqFNm9csV8OdDHGbAc+Aq4xxrznheMWLi4X9O4Na9favViTkqBbN2jaFL78MugDPjbWDssULQpt2mi4K+VLHge7iIwSkSoiUgO4CfhGRPp6XFlh5XJB376wfr2dFnnkiN1NukULO20yiAO+Zk0b7jExNtxXrXK6IqVCk85jD1Th4dC/P2zcCFOmwL59to3iFVfYhU9BGvAXXGCHZbLCffVqpytSKvR4NdhFJF5EOnvzmIVeRATccQds3gz/+x/s2GEHqlu1sq0LgkCGW1i8fj8TFm9m8fr91IgVvv0WoqJsuP/+u9MVKhVajDhw5RcXFyfLdUliwZw6Za/gn3kG9u6F1q3h//7PXskHoAy30G/aLyTsTCQlLYPoSBeNqpbi3QHN2fqHoXVr+yl9+y1cconT1SoV2IwxK0QkLrfn6VBMsClSBO66C/74A15+2e7LeuWVdnVrADZEj994gISdiSSnZSBAcloGCTsTid94gFq1bKBHRtruC2vXOl1tAHFnwMav4Lvn7Z/u0FrApnxLgz1YRUfDiBGwdavd4GPlSrvpR6dOAdWgZe2e4yTu3Myfz3Xm2M+fApCSlsG6PXZLwaxwDw+34b5unYcnDIVAdGfAu91hxu3w7TP2z3e75+9zCYWvgyowDfZgFxMDDzwA27bBuHH2qr1pU7u6NQBWA80eP5K9b90NQPHG1wEQHemiXqUS2c+pXdveUHW57MhSgcPdG4EYCDYvhN3L+X3ncUAgLQl2L7eP50WofB1UgWmwh4pixWDkSBvwTz9tu0o2bgw9ejgy9WTHjh0YY/jmy1lc2msEdR+bjysyipjMMfZWdcr/7fm1a9sr97Awe+W+viBdhzIDcd3u4xxOzsh/IAaKfat5eclRGryeRPLpzHtgacmwL493mTO/DqQlIeL27OugV/5BSYM91JQoAaNHw/bt8MQTdmpkw4bQq5ffBrEff/xxqlevDsD+/ftZ9uF4Xu19Kfe1q82rvS/l3QHNcYWZf72uTh0b7sbYK/cNG/J54n2rmfxTIhdPSmLNgcyWDPkJxACRcLQo936dysDGEcREZH6dImOgYv28HWDfakhLpuKLJ5i3ObN7aEG+DnrlH7Q02ENVyZIwZowN+EcfhfnzbaOx3r0LkJh5c/ToUYwxPPXUUwwfPhwRoXz58rjCDG3qVmB4m1q0qVshx1DPctFF8M039u/5DfeX5m9h8NwUbmsUztU1wu2D+QnEAJCUlMSl3e8CYHLPsoCByKJQOQ5qtfv3C3K6oq7YgLfWwP4kIa6Syz6vIF+HzCv/0yknSTnt4ZW/8isN9lBXujQ89ZQN+JEjYc4c2yq4Xz87N95Lpk6dynnnnQfAxo0bmTBhQoGPVbeuDXe324b7xo25v+bJJ5/k/uemcHe7WN68oRy5BmKAql/fhu/pU6nQczq0Hm3/7DcLwlx/f/LZrqgvuIaNqeW4uWEUFYq5Cv51yLzyf//305R/8YR9LAh/AyqURMTvb02aNBHlkP37RR54QCQ6WsTlEunfX+SPPwp8uNTUVImOjhZAOnbsKG6322ulrl0rUr68yPnni2zYcPbn3X///QLI6NGjRTLSRTbMF4l/3v6Zke61evxh3bp1cujQobw9ecN8kbHny5tdo6THReEiY0qIjD3/r8/b069D5vGTHykuX98c/ffjK0cAyyUPGavBXljt2ydy770iUVEi4eEid9whsn17vg4xb948AQSQpUuX+qTMNWtEypWz4b5x478/PmjQIAFk3LhxPjl/QIt/Tk49WkIAubKaywbvmJI2zL0hI13kretsmI8paf9867qg+2EZSjTYVd7s2SNy990ikZEiEREiQ4aI7NhxzpdkZGRIw4YNBZALL7xQ0tN9+x/9999tuFeqJLJp01+P9+7dWwB59dVXfXr+gLVhvsRVDhdAMh4v7psr6iD/DSjUaLCr/Nm5U+TOO224R0aKDBsmsnv3v562fPny7Kv0mTNn+q281atFypYVqVxZZPNmkU6dOgkgb775pt9qCDgZ6fLGbQ3kh0Fl9Iq6kMhrsGuvGPV3f/4JY8fCm2/aFUNDh9qbrhUr0qtXLz791K4eTUpKIiYmxq+lrV5t57ifOHGAtLSWfPLJs9xwww1+rSHP3Bl29si+1VCxgb1x+c+bn149z+921ouvzqMCQl57xWiwq5xlLXR6+23cERG8lJrK88BjEyYwfPhwx8qqVasnW7a8QdmyRfn552guuMCxUs4ua7bK7uV2FklkjJ2VktPMFqXyQZuAKc/ExsK0abwwYADvpqZyL7A/Jobhu3fDoUOOlFS1alW2bJnJ5MlbEYmmdWvbKifgZM7/HvDZIbp8mKTzv5XfabCrHB05cgRjDA9Nnszq++7DtWEDpnt3eP55G/qPPmp3d/IDEaFo0aLs2rWLpUuXMnBgMxYvtrsHtmoVgOG+bzUvxh9l+srTdKyVuVBK538rP9JgV//yxhtvUKZMGQC2bNnCf//7X7ve/733YM0a6NjRjsPHxtrVrYmJPqtFRAgLCyM5OZkVK1Zw2WWXAbZLwuLFcPKkXcS0bZvPSsi3matP8ODCVEY0j2BIXKR9MMhWwKrgpsGusqWmphIREcGQIUPo2rUrbrebC/45iF2vHnz8sb2T2bat3eQjNtaubj1+3Kv1uN1uwsLst+jatWtp3Ljx3z7eqJEN9xMnbLhv3+7V0xfIsmXL6Hn307SpV5aXu+ShJYBSvpCXqTPnegOqAt8C64G1wIjcXqPTHQPPnDlzsqcx/vLLL3l/4W+/iXTpYmfOli4t8swzIidOeFxPenp6dj1btmw553NXrLCnrl4932usvGr79u0CSKlSpXT+t/IJ/DWPHTgfaJz59+LAJqDeuV6jwR44MjIypF69egJI3bp1C77YaNkykU6d7LdU2bIizz8vcvJkgQ6VlpaWHeo7clkslWXFCpFSpURq1HAm3BMTE7Nr9mZbBaXOlNdg93goRkT2ishvmX8/kXnlXtnT4yrfW7ZsGS6Xi3Xr1vHFF1+wbt06XK4CTseLi4O5c+1GH02awEMPQc2aMH48pKTk61DTpk0DYN++fVStWjVPr2ncGBYtssP9rVvbPb/95fTp05QqVQqAtLQ0jDl790ql/CIv6Z/XN6AGsAMokcPHBgHLgeXVqlXz8c81dS5ut1u6d+8ugLhcLklOTvb+SX78UaRtW3sFX7GiyCuviKSk5Lm+jIyMAp122TKRkiVFYmNz7YzgNT169BBAEhMT/XNCVWjhryv2LMaYYsAM4B4R+dddNBGZLCJxIhJXrlw5b51W5dPWrVsJCwtj1qxZTJo0ifT0dKKjo71/opYtYeFC+O47O6NmxAi48EKYNAlOnTrnS40x2TdN8ysuzp728GE7FXLnzgIdJl+mTZvG4cOHKVmypO9PplQeeCXYjTER2FB/X0RmeuOYyvsefPDB7Fkuhw4dYujQob4/6VVX2W2RFi+GGjVg2DC7g/XkyZCW5pNTNm1qw/3QITsss2uXT06TrVSpUtm96JUKBB4Hu7EDitOA9SLykuclKV9Yt24dL774Ig899BAikj1P3S+MsU1evv8eFiyASpVg8GB7JT99Opw+7fVTNmtmT3XggL1y373b66dQKmB544r9cqAfcI0xJiHzraMXjqu8qG7duqSnp/Pcc885V4Qx0K4d/PQTzJsHZcvCgAF2y6R33oH0dK+ernlzDXdVOHljVswPImJEpIGINMp8m+eN4pT3GGMKPuPF24yBDh3g11/hiy+geHG49Va7Zd/770OG9zZLbtECvv4a9u+3wzIhFe457XeqfCPIvta68lSdna+/mY2B666D336DmTOhSBHo29duuv3xx3bTUy+47DL46ivYu9eOCO3Z45XDOuts+50GeOAEpX98rf/30I0YVzjbt/7hdGVnpcGucubP4DAGuneHhAT45BP7/k032YYwM2Z4JeBbtrThvmePvXLfu9cLdTsps4PkfXMO8cKPqdpB0pcyv9b7j5zAPHmMO784ztU1Iqh2Kg+7rDtEg13lbPNCErf8ihm9l3UH0/0THGFhcMMNtg/Nhx/am6rXX29XH82eDR7uHXD55TB/vh2OCfpw37ca0pIZ//NptiVm/uDTDpK+kfm13nHMTdUShh33FCP+1hjCDqx1urKz0mBXOUrc/Auln94PQM3Smd8m/goOl8tesa9dC+++a/vzdutm5zF++aVHAX/FFfbKfdcuOyyzb58X6/anig3YnmQ7Rz7TJso+ph0kfaNiA4iMoWnlcHbcW5yqJcMC/mutwa7+5dixY5S+7gkAUkYXJyo8c4m8v7+ZXS475r5+vd2q7+hR6Nz5rzuiBQz4K66wV+47d9or96AM91rtKBnbmDe6lqBUVJh2kPSlWu3s1zayKMHSrVO3xlN/c+zYsey+JymTOxJ1cGXgbO92+rSdFvnUU3Zv1pYt4cknoU0bOy6fT0uW2Nby1arZNVQVKvigZl/S/U79J0C+1rrnqRdluIX4jQdYu+c4F1cqQas65XGFhV6jp+PHj2cvi09JSSEqMiIgvpn/JS3NLmwaO9aOqVx5pe0L36pVvg/13Xc23GvUgG++CcJwd8jYsWPp1asXtWrVcrqUQkWD3Usy3EK/ab+QsDORlLQMoiNdNKpaincHNA+pcD8z1JOTk33TP8bbTp2CqVNtwO/da8dV/u//7FhLPsTHQ6dONty//RbKl/dJtSHFGEOrVq349ttvnS6lUNHNrL0kfuMBfl21js2vDcTtziA5LYOEnYnEbzzgdGleE5ShDnbe+7Bh8Mcf8PLLsG6dvXq/9lq7ujWPWrWy92S3bbM3VA8e9F3JoaJ27drEx8czYfFmFq/fT4bb/xeI6uw02HPx2sRX+eO1AaQf+auTVEpaBuv2eHcbOKecOHEiO9STkpKCJ9TPFB1tu0du3Qovvmjnw7dsacdYli3L0yFatbLt5LduzSXcg2wFoi9kuIVSzXsCMH7hJoZ/uJJ+037RcA8gGuxn4Xa7qVGjBrMmjaVkw3ZUHzkXkzm+HB3pol6lEg5X6LkTJ05QooT9PJKSkoiJiXG4Ig/FxMD999tL7+eesy0LmjWDLl1g5cpcX37NNTBnDmzZYu/HHjr0jyfoak/A/hZ7rEoLANKO7g3J32KDnQZ7Dv78809cLhd//vknCxYuouOw/yMm0oUBYjLH2FvVCe6B2JMnT4ZWqJ+paFG7g9O2bfD00/DDD3aRU/fusGrVOV/apo0N982bcwj3zBWI//vxKJ+sPVVoV3uu3XOcU24XRRtciwmPAELrt9hQoMH+D1OmTKFGjRqAHXtu17YN7w5ozqu9L+W+drV5tfelQX/j9OTJkxQvXjz77yEV6mcqXhxGj7YB/+STdtpLo0Z2devas68abNvW9ibbtMn+/fDhzA/sW038puPcOS+V3cczhx0K4WrPiyuVIDrSRdkOdxNevCwQOr/FhgoN9kwiQv369Rk0aBA33HADIpIdfq4wQ5u6FRjephZt6lYIqVAvWrToOZ+f4RYWr98f3DfJSpaExx+H7dvhscfs4qb69aF3b9iwIceXtGtnuxhs2PBXuB+IqEbrt5OoUcpw72VF7BMDfAWiL7SqU55GVUuF3G+xoUSnOwJ79uyhcmW7//a8efPo0KGDwxX5RlJSEsWKFQPs+HrW388mZKd6Hj5sb7K++qrdaLtPHxv8OczJ/vpr6NoV6tUTVq4sAxzF/XRFzOmUwFi05ZCstR3r9hynXgiv7Qg0Oo89j9577z369esHwJEjRyhdurTDFflGfkMdYPH6/fS+axQHv3mT6iPnAvbq7NXel9Kmbgis5Dl40Ab8xIl2Tny/fvaKvmbNvz3tq6+gQ4dUYB3b/oilxulfAm/RlioUdB57LkSEli1b0q9fP9q3b4+IhGyoJycn5zvUTypMercAACAASURBVJ8+zX8aVuXgN29SrNFfv8EE0k0yj4eJypWzs2e2boW774aPPoLateGOO+ywTaZPPx0AdCMiohE9byjNkXLt4eoHoU57DXUVkMKdLsAJBw4coELm2vEZM2bQo0cPhyvyneTk5Oxx9LyG+rJly2jWrBkANQa8jJS9MPtjgXKTzKvDRBUqwEsvwYMPwrhx8PrrtifNgAF82bAh06dPZ8KECVx4YRjdutn1TwsXQoheB6gQ4JUrdmNMe2PMRmPMFmPMw944pq/MmDEjO9QPHDhQaEL9+PHjeQr1wYMH06xZM4oUKUJK6ikua948IG+SxW88QMLORJLTMhDwzlzq88+HV16xK1kHDkSmTaPt0KHMqlyZ4T170qEDzJoFv/9uwz0x0WufjlLeJSIevQEu4A+gJhAJrALqnes1TZo0ESdce+21Asjll18ubrfbkRr8JSkpSQAB5Pjx47k+/9ixY9nPHzduXPbj6RluWbRun0xYtEkWrdsn6RmB8XV7ZdEmKXVlXwGk0sA3pPrIuVJj5FyZsGiTV46fmpoq1UBeB5HwcJEiRURGjBDZu1fmzBGJiBCJixM5etQrp1MqT4Dlkodc9sYVezNgi4hsFZE04COgqxeO6zVHjhzBGMOCBQt4//33+eGHHzAFaPMaLFJSUrKv1I8dO5Y9vfFsvvzyy+y2Alu2bGHkyJHZHwvEqZ4HDx5kRNvaJH7/HkUbXEvEeXZGkzeHiaKiotgB3Jaaaie033yzvclasyad4x9g5vREVq3SK3cVmLwR7JWBnWe8vyvzsb8xxgwyxiw3xiw/6McuS19++SVlypQBYPfu3fTp08dv53ZCSkpK9oKjY8eOZa8uzYmI0LZtWzp37kz9+vVxu91ccMEF/iq1QMaOHUv5zPaLnZ/+lGpd7/X6MJHb7aZVq1Zs3ryZIkWKQGwsTJtmJ7XfcAOMH0/nwZWZcd1bJCQI//kPHDvm8WmV8p68XNaf6w24AZh6xvv9gFfP9Rp/DcVkZGQIIA0bNgz5oRcRkZSUlOzhlMTExHM+d/fu3dnPfe+99/xUYcEdOHAgu94BAwaIiIPDRBs2iPTpI2KMzI7qJRFhp6V5k9OSy5dcKY/hx6GYXUDVM96vAuzxwnE9FhYWxqlTp0hISAjpoReA1NTU7M6MiYmJ2UMrOZkyZUr2gqwDBw5w8803+6XGgnrmmWeyr9I3b97M1KlTAQeHierUgfffhzVr6NLV8Km7JytWCO0v2cnxHaE3LiMinDp1yukyVH7kJf3P9YadMrkViOWvm6cXn+s1Tt08DVVnXqkfPcfdvPT0dKlataoA0qNHDz9WWDBnXqX379/f6XLObvVqmdVinISTJpe5fpFjo58TOXbM6ao85na75eGHHxZAWrdu7XQ5SvJ+xe5xsNtz0RHYhJ0dMzq352uwe09eQ33t2rXZz1u8eLEfKyyYZ555JrvezZs3O11Onsx4cauEm9PSkh/keOlqIs88I3LihNNlFciZX/8uXbrI6dOnnS5JiZ+DPb9vvgj2rPHWVwJsWp4vpaamZv/nO3LkyFmf98gjj2Q/LykpyY8V5t/Bgweza73tttucLiffPvtMxOVyy+Wl18hxiomULSvy/PMiJ086XVrOMtJFNswXiX9OZMN8mfjqhOyvf6tWrSQ1NdXpClUmt9tduII9PcMtvSf/JHUfmy81Rs6Vuo/Nl96Tfwr5cJ81a9Y5Q/3Mq/kHHnjAz9Xl37PPPptd76ZN3pmP7gQb7iJXNDwmx9t0s//NypcX+e9/RZKTnS7vLxnpIm9dJzL2fHm7W3T2175Ro4ZyMlB/EBVSf/zxR9a/j99unjouaxVi0qn0kN2XNCfdunVDJOceN0uWLMm+mbpq1SpeeOEFf5eXZ4cOHcIYw6hRo7jtttsQEWrl0GkxWPTsadvO/LSmBJ3SZnFywVLbJvj++22DsQkTIDXV6TKzNw6Z8tMRbv08hdhShqOPVmDlR+Nybees/Gfs2LH5noYcEsG+ds9xjm5dzY7nr+PYz58CgdWsyt/69OnD1VdfTfny5UlPT6dBgwZOl3RWz44bR7ly5QB4a95Spk6b7nBF3nH99fDBB7B0KXR6+jJOfr4IvvvOzqgZMQIuvBAmTbJdJZ2ybzWkJXN9vUgOP1ScrSOKU8qVWug2DglUKSkpGGN49NFHGTp0aNb9zDwJ+mAXET546k72vT8SU6QoJZtfDwROsyp/ylph++GHHzJx4kT279+PyxWY3QcPHz6MMYZHRo2iRIO21Bg5l+d+OhZSmyL36mVnRf7wA3TqBElNroL4eFi82C56GjbM9oCfPBnS0vxfYMUGEBlD6WjDedGZU0UL4cYhgWjx4sXZCw0TEhKYNGlSvl4f1MG+fft2wsLC+OW7RTTrN4qLHvqMMFd4QDWr8pdVq1Zlr7DdsWMHw4YNc7iic8saGqp55xSi6rXm6A8fhOQQ2o03/iPck7C7Zi9ZAgsWQOXKMHiwvZKfPh1On/ZfcbXa2Y1CIosCxv5ZOc4+rvzHnQEbv4Lvnkc2zKdD+/a0bduW6tWrk56eTsOGDfN/zLwMxHv7zRs3T+Pj47Nv9hw6dChgm1X5y/79++WBBx4IqhW24xdskMjyNQWQmLpXeb2RVyD54AORsDCRVq3+MUHG7RaZP992FAORCy4QeestEX9NL8yeFfO8/TMj3T/nVdYZN7B33ls8O9Pemj49x6cT6rNiDh48KJMmTfL4OMoZS5cuzf4mLtttlFQfOVeqZ85oWrRun9Pl+cT779twb91a5F+zTt1ukS++EGnUyP63rF3bviBdgzakbZgvMvZ8kTElsv8/HHykgn08B3kN9qAdiilbtixDhw51ugyVTyJC69atadmyJdHR0fR6LZ5y9a8KuH7vvtCnj92/47vv4LrrIDn5jA8aYx/87TeYOROiomxHyfr14eOPwe12rG7lQ5k3sAG23l0MGVOCshGe38AO2mBXwSchIYGwsDDi4+N5++23SU5O5oMhV/Fq70u5r11tXu19afBvlJ2Lm2+Gt9+Gb7+FLl3+Ee5gA757d1i5Ej75xL5/003QsCHMmKEBH2oyb2ADxJbOjGMv3MDWYFc+JyJ0796dSy+9FLBb9N1yyy1AYPZ797W+fW24f/MNdO0KKSk5PCkszLYIXr0aPvzQ3lS9/npo3Bhmz4Z8TH1TAcxHN7A12JVPbdiwgbCwMD7//HMmTpyIiORpi75Q168fvPmmnfnYpctZwh3A5bJX7GvXwnvv2Uv8bt0gLg6+/FIDPtiFuaDfLOg5HVqPtn/2m+XxJuka7Mpn+vfvT926dQE7xz7Qp2D626232hmOixef48o9i8tlx3HWrYO33oKjR6FzZ2jRAr7+WgM+mIW5oE57uPpB+6eHoQ4a7MoHtm3bhjGGt956i3Hjxp217YGC226zmzMtWmSH1nPtNBAebn8ibNwIU6bAvn3Qvj1ccYU9iAa8QoNdedn48eOpWbMmAPv37//b/qkqZ/37w9Sp9sK7W7c8tpGJiIA77oDNm+F//4MdO6BdO7j6aru6VRVqGuzKq6pXr87TTz+NiGTveqRyd/vtf4V7jx756BEWGQlDhsCWLXaz7T/+gNat7erW77/3ac0qcBlx4Fe3uLg4Wb58ud/Pq1SgmzoVBg6Ejh3tdPYiRfJ5gNRUeOMNePZZ2L8f2raF//s/uOwyn9Sr/MsYs0JE4nJ7nl6xOyzDLSxev58JizezeP3+gGiAFYg1FRZ33GF7gs2bZ6/c8938MSrKdo/cuhVefBFWrYKWLaFDB/j1V5/UrAKPXrE7KMMt9Jv2Cwk7E0lJyyA6c+Wlk4t0ArGmwmjyZNsbrFMnuy4p31fuWU6ehNdeg+efhyNH7EyaJ5+08+FV0PHLFbsx5gVjzAZjzGpjzCxjTClPjlfYxG88wMo/D3N4828IBER3w/iNB1i54yh7l3zC6ROHA6KmwmjQIHj9dTtV/YYbPOjqW6wYjBwJ27fD00/Djz9CkyZ2Cs6qVd4sWQUQT4diFgKXiEgD7GbWozwvqfD4bPaXbBjbmQMfjc5+zOkNQhb+uJwNT3ck8bs3bTvRAKipsBo82O7FMWeOh+EOULw4jB4N27bZK/Zvv4VGjeyB16zxWs0qMHgU7CKyQETSM9/9GajieUmhLz09nVq1avH6qAGElyhHtQc+z/6YUxuEiAg9evTguYGdAag64iPCS5Z3tCYFQ4fakZQvvrAbd3i8H0fJkvD44zbgH33UTsNp0AB694b1671Ss3KeN2+e3g7MP9sHjTGDjDHLjTHLDx486MXTBpcffviBiIgItmzZwuezv+CGF7+gaHQRR7sb/v7774SFhTFr1ixefuUVek/+iWIlShaKjovB4M477UzG2bPtxh1e2WypdGl46ikb8A8/bH8tuPhi28hm0yYvnEA5Kre+vsAiYE0Ob13PeM5oYBaZN2Nze/NGP/ZgNXHiRLnoooskLS1NRMTRDULcbrd07949uw90YmKi4zWps5swwbZq795dJPPbx3sOHBB58EGR6GjbNP7WW0W2bPHySZSn8NdGG8CtwE9ATF5fU5iDPZBcc801Asirr77qdCnqH7755hsZO3bsvx5/5RX7v7ZHDx+Eu4jIvn0i994rEhUl4nKJDBggsm2bD06kCiKvwe7RdEdjTHvgJeBqEcnz+IpOdwwMe/fupWjRopQooePngURECAsLo27duqxbt+5fH3/5Zbj3XujZ03b0jYjwQRF798K4cXZqjohdGjt6NFSt6oOTqbzy1wKliUBxYKExJsEY87qHx1N+dP7552uoB6DBgwcDcLaLn3vugZdesvPb+/Tx0f7X558Pr7xiWxTccYdtQ3nhhXDXXbB7tw9OqLzJ01kxF4pIVRFplPk2xFuFKVWYZK32fWFuAlOmTGH43XcTExNz1uffey/897/w2We2m296+lmf6pkqVeycy82bbSvKN96ACy6wq1v37vXRSZWntKWAUg7LWu07/MOVjLrpGgAOXnxTrq0c7rvPdg349FMfhztA9eo21Ddtsid77TWoWRPuvx8O6OK1QKPBrpTD4jceIGFnIslpGWQkHaVcj8dYtetYnlb73n8/vPCC3R61b18fhztAbKxtIL9hg51Y//LL9rGRI+HQIR+fXOWVBrtSDlu75zgpaXaVb7UHZxNTq3m+Vvs+8IBtBfPxx3bLPZ+HO9jx9rfftjs6de9uf7rExtobrEeO+KEAdS4a7Eo57OJKJYiOtNuhmcxt0fK72vfBB+G55+Cjj+wGS34Jd4A6dexerGvW2F7Dzz5rA37MGEhM9FMR6p802JVyWKs65WlUtRQxkS6PVvs+9JDN1Q8+sOGekeGbenNUr579lWHVKruT0//9H9SoYVe3Htc+Q/6mbXuVCgAZbiF+4wHW7TlOvUolaFWnfIHbJD/7LDzyiL3H+fbbdh9sv0tIgCeesH0QSpe240XDh9tmZKrA8jqPXYNdqRD0zDN2uLtfP3jzTYfCHWDFCjss8+WXULas/bXizjuhaFGHCgpuuoOSUoXYI4/YUZB337WLRv06LHOmJk1g7lz4+WeIi7PBXrOmXWGVkuJQUaFPg10FFN2Wz3sefdQOdb/zDgwY4GC4AzRvDvPn240+6te38zRr1oQJE/Kxc7fKq3CnC1Aqi27L532PPQZutx3uNsZulu3YsAzY/VcXLYIlS+wQzYgRdq7mI4/Ynz4F3gNQnUmv2FXA+CphO7PHDmbDuO4Bs1VgKBgzxr699RYMHGiD3nFXXWV3cfrmGzs9ctgwqFXLrm71SsP5wk2DXQWEDh060LlJTZL/XE3JFjdkP67b8nnHE0/YjZPefDOAwh2gdWt79b5gAVSuDEOG2Lnx06f7qLtZ4aDBrgJC6dKlGTLySeo+Np+SLW/Mfly35fOeJ56w4+7Tp9v9VAMm3I2xc9+XLv1r9syAAVC3rr1B4LfVVqFDg10FhA8++ICJzzzmlYU6KmfG2Jupo0fbsfYhQwIo3MEW2LEj/Pqr3eS1RAm70qpePXj/fYfv/gYXnceuAoo3F+qonInYK/dnnrFX7pMmQVggXuKJ2AVOY8bA6tVw0UX2144bbgjQgn1P57EHGZ3mZ7nCDG3qVmB4m1q0qVtBQ90HjIGnn4ZRo+y9ymHDAuzKPYsx0K0brFxpexO7XHDTTdCwod1lJCCLDgw63TEA6DQ/5W/GwNixNhufe86+/9pr9s+AExYG118PPXrYgH/iCft+w4bw5JPQpUuAFu4cvWIPAPEbD/DLilVsffsh0lOO6zQ/5RfG2L4yDz0E//ufvXJ3YGQ278LC4MYbbSfJd9+F5GR7RR8XZ2+6BnTx/uWVYDfGPGCMEWNMWW8cr7BZu+c4W6cO59SO3wmLiAZ0ml+BuTNg41fw3fP2T7fecDsXY+ye1Q8+aMN9+PAgyEeXy+4qsm6dnb959Ch07gwtWsBXXwXBJ+B7Hg/FGGOqAu2AHZ6XUzhdXKkEF42aTcrpv8YMdZpf/mS4hfj1e6m14BYqJa3FlZ6CiYyBynHQbxaEObncMrAZY4djROxWe8bYlf4BP7oRHm73Yb35Zjst8qmnoEMHuOwyO/2nTZsg+CR8wxtX7OOBhwD9MVlAreqU59JqpXWaXwFl3aP47OM3KZO4mvD0ZOLeOMFlk/aTvG0ZbF7odIkBzxi7sv/++2HiRLvSP2gufCMi7Lz3TZvg9ddh5047L/7qqyE+3unqnCEiBX4DugCvZP59O1D2HM8dBCwHllerVk3U36VnuGXRun0yYdEmWbRun6RnuJ0uKWgsWrdPKnV7SDq1jpO0x0qIjCkhkzpGCfZiQy67uJokJSU5XWZQcLtF7rtPBERGjLDvB53UVJGJE0UqVbKfSKtWIt9/73RVXgEsl7xkc65PgEXAmhzeugK/ACUlD8F+5luTJk389GVQhcHLCzdKWFRxAWTng2VExpTIfnv22uLZAb9kyRKnSw0KbrfIPffYdLjnniANdxGR5GSR8eNFKlSwn0y7diJLlzpdlUfyGuy5DsWISFsRueSfb8BWIBZYZYzZDlQBfjPGVPTsdwil8ueSyiWp8+AnxI78gu1FapMkRcgQQ3p4DA/3aYVkpPPWW2/RqFEjp0sNCsbYdukjRsDLL9vhmaAZljlTdDTccw9s3WpvHiQk2O6SHTvCsmVOV+dTXlt5mhnucSJyKLfn6spT5U1nrgM4lXaaayNX07rkfnp26oCr9rV647SARGwuTpgA9933143VoJWUZG8gvPACHD5sZ9I8+SQ0bux0ZXnm963xNNiVk7QVgW+I2Cv3V1+125Y+/3yQhzvAiRP2E3rxRTtVsls3u+ipYUOnK8uV7nmqlPIKEbj7bnux++CDf61UDXrHjsErr9hxp2PHoGdPG/CXXOJ0ZWelvWKUUl6RNa/9zjvtKMaoUUE65v5PJUvaJvXbttmuaAsWQIMG0Ls3bNjgdHUe0WBXSuXKGHvFPnSovWJ/5JEQCXeA0qXt4qZt2+Dhh2HOHLj4YujXDzZvdrq6AtFgV0rlSVa4Dxli2xCMHh1C4Q5QpoztZbxtm50KNHOm3eyjf387syaIaLArpfIsLMx2gRw0yDYQe/TREAt3gHLl7F3irVvtneOPPoLateGOO2D7dqeryxMNdqVUvoSF2YZhAwfaC9zHHw/BcAeoUAH++18b8MOGwXvv2Q23hwyxbQsCmAa7UirfwsJsW5Y77rCbdowZE6LhDnD++Xb2zJYt9qfZ9Olw4YVw112we7fT1eVIg10pVSBhYXYHpgED7L3HJ590uiIfq1LF7iO4ebPtKvnGG3DBBXYV1759Tlf3NxrsSqkCCwuDyZPh9tttsId8uANUr25DfdMm2zJ44kSoWdPecD0QGJvjaLArpTwSFgZTptjJI088YVuhFwqxsTBtmp3z3quXbawTGwsjR8KhXBfg+5QGu1LKY2FhMHWqHaEYM8YOzRQaF14Ib70F69dD9+52FVdsrJ0PeuSIIyVpsCulvCIr3G+5xc6UGTv27x/fvn07l19+Oenp6c4U6Gu1a9uZM2vWQKdOdspQjRr2J11iol9L0WBXSnmNy2UnjfTrZ+e4P/OMfXzDhg3ExsaydOlSXK4Q77ZZr56d+756NVx7rR2bio21v8Yc988+xhrsSimvcrnsHtN9+9rRiOHD91C3bl0A3G43JiQ6iOVB/frw2WewcqXdpu/xx23AP/us7TDpQxrsSgWQkydPOl2CV7hcdtj52msPMnFiJVyuRwpXqJ+pUSP4/HNYvtxutP3IIzbgn3/e9oj3AQ12pQLE119/TfHixTkQIFPmPBUfv5gFCyoSFTWTjIyxPP98IQz1MzVpAnPnwi+/QNOmdvZMzZq2bXBysldPpcGulIMy3MLi9ft58sMltG/fngYNGlK+fHmny/LYnDlzaNu2LRdfXJcTJ3rQu7dtnPjCC05XFgCaNYP58+HHH22b4PvvtwudJkyA1FSvnEKDXSmHZG3pd9d7y3iiz9UA1Bv2PzLcwb02/6OPPqJLly5ceeWVrFmzhvBweOcduOkmeOghu3GRwu6/unAhfPcd1KljG45dcIHtsnbqlEeH1mBXyiHxGw+QsDORw+t/AqDqfTNYtesY8RuDdyhm6tSp9O7dm65du7JkyZLsx8PD4d134cYb7S5M//2vg0UGmquugvh4+OYbOzRz11222djkyZCWVqBDehzsxpjhxpiNxpi1xpjnPT2eUoXF2j3HSUnLIPrC5lR7aA5hEUVISctg3R7/TInztvHjxzNw4EBuueUWPv/88399PDzcTvPu1cvunzp+vANFBrLWrWHJEruTU+XKMHiwvZKfPh1On87XocI9qcMY0xroCjQQkVPGmOAfHFTKTy6uVILoSBfJZ1yURUe6qFephHNF5UPWBuJr9xxnxazJvPPaiwwfPpwJEyac9TXh4fD++7YT5H332cfuvddPBQcDY6BdO2jbFr7+Gh57zHZZy+qPnEceBTswFBgnIqcARCR4f4dUys9a1SlPo6qlSNiZaK/cI100qlqKVnUC//oo6/5Aws5Eds1/gxPLP6dex9sY//Irub42K9zdbhvuxtgGieoMxkD79vCf/8CXX9pQv/XWPL/c02CvDVxpjBkLpAIPiMgyD4+pVKHgCjO8O6A58RsPsG7PcepVKkGrOuVxhQX+tMCs+wN7vvuQE8s/p1Sr/kiTXsRvPECbuhVyfX1EBHz4ob2heu+9th3B3Xf7ofBgYwx07mxbFMyebXvR5EGuwW6MWQRUzOFDozNfXxpoATQFPjHG1BT5d8t9Y8wgYBBAtWrV8lScUqHOFWZoU7dCnsIwkGTdHwgvXobz/nMXxRu1z74/kNfPJSLCrry/8UY7ISQszN43VDkwBrp1y/PTcw12EWl79nOZocDMzCD/1RjjBsoCB3M4zmRgMkBcXFxwz+dSqpDLuj8gl7TJfqwg9weywr1XLxg+3ObXsGHerrbw8XRWzOfANQDGmNpAJOBsI2KllM9l3R+IiXRhgBgP7g9ERsInn0DXrvaKfdIk79db2Hg6xj4dmG6MWQOkAbfmNAyjlAot3r4/kBXuN9xgr9iNgaFDvVx0IeJRsItIGtDXS7UopYKIt+8PREbCp5/C9dfDnXfacB8yxCuHLnQ8vWJXKiCdOcf64iCabVLYnRnuQ4facB882Omqgo8Guwo5WXOsV+44Suppd/b88HcHNNdwDwJFitg25j172it2Y2DQIKerCi7aK0aFnKw51n9+9ASp+7aQnJZBws7EoO7BUtgUKQIzZkDHjvaKfcoUpysKLhrsKuRkzbEufXV/IivUBAjqHiyFVVa4d+hgr9inTnW6ouChwa5CTtYc64iyVTHGfosHUw8W9ZeoKJg5066uHzgQpk1zuqLgoMGuQo4351gr50VFwaxZtm3KwIF2P1V1bnrzVIWcYO7BonIWFWW3De3a1TY7NAZuu83pqgKXBrsKScHag0Wd3ZnhfvvtNtzz0fCwUNGhGKVU0IiOtk0O27SB/v3trkzq3zTYlVJBJSvcr7nGXrG/957TFQUeDXalVNCJiYEvvrC7yd16q924Q/1Fg10pFZRiYmDOHLj6arjlFvjgA6crChwa7EqpoJUV7lddBf362V2ZlAa7UirIFS0Kc+fClVdC3752447CToNdKRX0iha1ez5fcQXcfDN8/LHTFTlLg10pFRKywv3yy224f/KJ0xU5R4NdKRUyihWDefOgZUvo08f2di+MNNiVUiElK9wvuwx697a93QsbDXalVMjJCvcWLeCmm2z738LEo2A3xjQyxvxsjEkwxiw3xjTzVmFKKeWJ4sVh/nxo3tyG+8yZTlfkP55esT8PPCkijYDHM99XSqmAkBXuTZvCjTfa9r+FgafBLkDW7gUlgT0eHk8ppbyqRAn46isb7r162Q6Roc6ISMFfbExd4GvAYH9ItBSRP8/y3EFA1pa0lwBrCnxi/ykLHHK6iDzQOr0nGGoErdPbgqXOOiJSPLcn5RrsxphFQMUcPjQaaAN8JyIzjDG9gEEi0jbXkxqzXETicnue07RO7wqGOoOhRtA6vS3U6sx1o41zBbUx5h1gROa7nwK63axSSjnM0zH2PcDVmX+/Btjs4fGUUkp5yNOt8QYCrxhjwoFU/hpDz81kD8/rL1qndwVDncFQI2id3hZSdXp081QppVTg0ZWnSikVYjTYlVIqxDgW7MHUjsAYM9wYs9EYs9YYE7Cra40xDxhjxBhT1ulacmKMecEYs8EYs9oYM8sYU8rpms5kjGmf+e+8xRjzsNP15MQYU9UY860xZn3m9+OI3F/lDGOMyxiz0hgz1+lazsYYU8oY81nm9+V6Y8xlTteUE2PMvZn/3muMMR8aY6LO9Xwnr9iDoh2BMaY10BVoICIXAy86XFKOjDFVgXbADqdrOYeFwCUi0gDYBIxyuJ5sxhgX8BrQAagHyAQjXwAAA0VJREFU9DbG1HO2qhylA/eLSF2gBTAsQOsEOxV6vdNF5OIV4CsRuQhoSADWa4ypDNwNxInIJYALuOlcr3Ey2IOlHcFQYJyInAIQkQMO13M244GHsF/XgCQiC0QkPfPdn4EqTtbzD82ALSKyVUTSgI+wP9ADiojsFZHfMv9+AhtElZ2t6t+MMVWATgTw2hZjTAngKmAagIikiUiis1WdVTgQnTkDMYZc8tLJYL8HeMEYsxN7FRwwV2//UBu40hjzizHmO2NMU6cL+idjTBdgt4iscrqWfLgdmO90EWeoDOw84/1dBGBgnskYUwO4FPjF2Upy9DL2QsPtdCHnUBM4CLyZOWQ01RhT1Omi/klEdmMzcgewFzgmIgvO9RpP57GfUx7aEdx7RjuCaUCu7Qh8IZc6w4HS2F97mwKfGGNqip/nieZS4yPAtf6s52zOVaeIzM58zmjskML7/qwtFyaHxwL2tx9jTDFgBnCPiBx3up4zGWM6AwdEZIUxppXT9ZxDONAYGC4ivxhjXgEeBh5ztqy/M8aUxv72GAskAp8aY/qKyHtne41Pgz1Y2hHkUudQYGZmkP9qjHFjGwYd9Fd9cPYajTH1sf/gq4wxYIc3fjPGNBORfX4sETj31xLAGHMr0Blo4+8fjrnYBVQ94/0qBOjwoDEmAhvq74tIIHYZvxzoYozpCEQBJYwx74lIX4fr+qddwC4RyfqN5zNssAeatsA2ETkIYIyZCbQEzhrsTg7FBEs7gs+x9WGMqQ1EEkBd4ETkdxEpLyI1RKQG9pu1sROhnhtjTHtgJNBFRJKdrucflgG1jDGxxphI7M2pLxyu6V+M/ek9DVgvIi85XU9ORGSUiFTJ/H68CfgmAEOdzP8jO40xdTIfagOsc7Cks9kBtDDGxGT++7chl5u8Pr1iz0VB2xH423RgujFmDZAG3BpgV5rBZCJQBFiY+dvFzyIyxNmSLBFJN8bchW1D7QKmi8hah8vKyeVAP+B3Y0xC5mOPiMg8B2sKZsOB9zN/mG8F+jtcz79kDhN9BvyGHcJcSS6tBbSlgFJKhRhdeaqUUiFGg10ppUKMBrtSSoUYDXallAoxGuxKKRViNNiVUirEaLArpVSI+X9Ov9Ll8MvdFQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "lr = LogisticRegression(solver='lbfgs')\n",
    "lr.fit(Xtest.numpy(), ytest.numpy())\n",
    "beta_train = slope.numpy().flatten()\n",
    "beta_test = lr.coef_.flatten()\n",
    "b_train = intercept[0, 0].numpy()\n",
    "b_test = lr.intercept_[0]\n",
    "hyperplane = lambda x, beta, b: - (b + beta[0] * x) / beta[1]\n",
    "\n",
    "Xtrain_np = Xtrain.numpy()\n",
    "Xtrain_grad_np = Xtrain_grad.numpy()\n",
    "ytrain_np = ytrain.numpy().astype(np.bool)\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(Xtrain_np[ytrain_np, 0], Xtrain_np[ytrain_np, 1], s=25)\n",
    "plt.scatter(Xtrain_np[~ytrain_np, 0], Xtrain_np[~ytrain_np, 1], s=25)\n",
    "\n",
    "for i in range(m):\n",
    "    plt.arrow(Xtrain_np[i, 0], Xtrain_np[i, 1],\n",
    "              Xtrain_grad_np[i, 0], Xtrain_grad_np[i, 1], color='black')\n",
    "\n",
    "plt.xlim(-8, 8)\n",
    "plt.ylim(-8, 8)\n",
    "\n",
    "plt.plot(np.linspace(-8, 8, 100),\n",
    "         [hyperplane(x, beta_train, b_train)\n",
    "          for x in np.linspace(-8, 8, 100)], color='red', label='train')\n",
    "plt.plot(np.linspace(-8, 8, 100),\n",
    "         [hyperplane(x, beta_test, b_test)\n",
    "         for x in np.linspace(-8, 8, 100)], color='blue', label='test')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
